#################################################################
##DO SELF-REPORTING REGIMES MATTER? EVIDENCE FROM THE CAT########
## Cosette D. Creamer and Beth A. Simmons #######################
## February 2019 ################################################
### Replication Code for Appendix A #############################
#################################################################
##### NOTE: this code uses output generated by running the ######
#######scripts: "CAT_CIRI MSM Analysis.R" and ###################
#########  "CAT_LHRPS MSM Analysis.R" ###########################
## Clear environment
rm(list=ls())
sink("CAT_AppendixA")
### Libraries
library(data.table)
library(MASS)
library(car)
library(SparseM)
library(stargazer)

##############################################
### Data loading and pre-processing
##############################################


## Load CAT data (full)
catfull <- read.csv("CATReporting_Final.csv", header=TRUE)

## Subset data to only include observations for CAT ratifiers, for observations > t+1 (t = year of ratification)
cat <- subset(catfull, catfull$CATob==1)
cat$CATregiondensityl1 <- 100*cat$CATregiondensityl1
cat$HRTpercentl1 <- 100*cat$hrtpropl1
cat$yr <-as.factor(cat$year)
cat$lngdpl1 <- log(cat$gdpl1)
cat$lnpop <- log(cat$poptot)
## Remove observations with missingness for covariates used for each Weighting Model

cat_wm1 <- subset(cat, !is.na(cat$lngdpl1) & !is.na(cat$lnpop) &  !is.na(cat$farissl1) &
                !is.na(cat$polityl1) & !is.na(cat$CATreview)  & !is.na(cat$nhril1) 
              & !is.na(cat$CATart22l1)  & !is.na(cat$CATregiondensityl1) & !is.na(cat$nevdem) &
                !is.na(cat$regionall1) & !is.na(cat$transdeml1) & 
                !is.na(cat$transition2l1) &!is.na(cat$HRTpercentl1))

cat_wm2 <- subset(cat, !is.na(cat$lngdpl1) & !is.na(cat$lnpop) & !is.na(cat$tortl1)  & 
                !is.na(cat$polityl1) & !is.na(cat$CATreview)  & !is.na(cat$nhril1) 
              & !is.na(cat$CATart22l1)  & !is.na(cat$CATregiondensityl1) & !is.na(cat$nevdem) &
                !is.na(cat$regionall1) & !is.na(cat$transdeml1) & 
                !is.na(cat$transition2l1)  &!is.na(cat$HRTpercentl1))

cat_wm3 <-subset(cat, !is.na(cat$lngdpl1) & !is.na(cat$lnpop) &  !is.na(cat$farissl1) &
                   !is.na(cat$polityl1) & !is.na(cat$CATreview)  & !is.na(cat$nhril1) 
                 & !is.na(cat$CATart22l1)  & !is.na(cat$CATregiondensityl1) & !is.na(cat$nevdem) &
                   !is.na(cat$regionall1) & !is.na(cat$transdeml1) & 
                   !is.na(cat$transition2l1) &!is.na(cat$HRTpercentl1) & !is.na(cat$transparencyindex_l1))

cat_wm4 <- subset(cat, !is.na(cat$lngdpl1) & !is.na(cat$lnpop) & !is.na(cat$tortl1)  & 
                    !is.na(cat$polityl1) & !is.na(cat$CATreview)  & !is.na(cat$nhril1) 
                  & !is.na(cat$CATart22l1)  & !is.na(cat$CATregiondensityl1) & !is.na(cat$nevdem) &
                    !is.na(cat$regionall1) & !is.na(cat$transdeml1) & 
                    !is.na(cat$transition2l1) & !is.na(cat$transparencyindex_l1) &!is.na(cat$HRTpercentl1))

##############################################
### Weighting Models 1-4 (Appendix Table A1)
##############################################

wm1 <- glm(CATreview ~ lngdpl1 + lnpop +nhril1 +  polityl1+ nevdem + transdeml1+ transition2l1+  
             farissl1 + HRTpercentl1 + CATart22l1 +CATregiondensityl1+  regionall1+
             CATnumprevrevs+CATrevprevyr + CATyrssincereview + yr,family="binomial", data=cat_wm1)

wm2 <- glm(CATreview ~ lngdpl1 + lnpop +nhril1 +  polityl1+ nevdem + transdeml1+ transition2l1+ 
             tortl1 + HRTpercentl1 + CATart22l1 + CATregiondensityl1+ regionall1+
             CATnumprevrevs+ CATrevprevyr + CATyrssincereview + yr,family="binomial", data=cat_wm2)

wm3 <- glm(CATreview ~ lngdpl1 + lnpop +nhril1 + transparencyindex_l1 + polityl1+ nevdem + transdeml1+ transition2l1+  
             farissl1 + HRTpercentl1 + CATart22l1 +CATregiondensityl1+  regionall1+
             CATnumprevrevs+CATrevprevyr + CATyrssincereview + yr,family="binomial", data=cat_wm3)

wm4 <- glm(CATreview ~ lngdpl1 + lnpop +nhril1 + transparencyindex_l1 + polityl1+ nevdem + transdeml1+ transition2l1+ 
             tortl1 + HRTpercentl1 + CATart22l1 + CATregiondensityl1+ regionall1+
             CATnumprevrevs+ CATrevprevyr + CATyrssincereview + yr,family="binomial", data=cat_wm4)

stargazer(wm1, wm2, wm3, wm4)


##########################################################################
## Table A2: Effect of CmAT Review and Review History (Appendix Table A2) ##
##########################################################################

### Weighting Model 1 ###
# Denominator model (YEAR + treatment(review) history variables + time-varying covariates)
denom.p <- predict(wm1, type = "response")
# Numerator model (YEAR + treatment(review) history variables)
numer.fit <-  glm(CATreview ~  CATrevprevyr + CATyrssincereview + CATnumprevrevs + 
                    nevdem + yr, family="binomial", data=cat_wm1)
numer.p <- predict(numer.fit, type = "response")
## Calculate weights
## Calculate weights
d.pr.tr <- ifelse(cat_wm1$CATreview == 1, denom.p, 1-denom.p)
n.pr.tr <- ifelse(cat_wm1$CATreview == 1, numer.p, 1-numer.p)
cat_wm1$w <- n.pr.tr/d.pr.tr
summary(cat_wm1$w)
sd(cat_wm1$w)
wm1dt = data.table(cat_wm1)
wm1dt[ , swrevf:=NA_real_]
wm1dt[ intersect(order(year), which(!is.na(w))), swrevf := cumprod(w), by=cowcode]
cat_wm1 <- as.data.frame(wm1dt)
summary(cat_wm1$swrevf)
sd(cat_wm1$swrevf)

#### Outcome MSM with LHRPS (t+2), Weighting Model 1#######

msm1 <- lm(farissf2 ~ CATreview + CATnumprevrevs+
            CATyrssincereview + CATrevprevyr+  CATlintt, weights = swrevf, 
          data = cat_wm1)

############################################################
### Weighting Model 3 ###
# Denominator model (YEAR + treatment(review) history variables + time-varying covariates)
denom.p <- predict(wm3, type = "response")
# Numerator model (YEAR + treatment(review) history variables)
numer.fit <-  glm(CATreview ~  CATrevprevyr + CATyrssincereview + CATnumprevrevs + 
                    nevdem + yr, family="binomial", data=cat_wm3)
numer.p <- predict(numer.fit, type = "response")
## Calculate weights
## Calculate weights
d.pr.tr <- ifelse(cat_wm3$CATreview == 1, denom.p, 1-denom.p)
n.pr.tr <- ifelse(cat_wm3$CATreview == 1, numer.p, 1-numer.p)
cat_wm3$w <- n.pr.tr/d.pr.tr
summary(cat_wm3$w)
sd(cat_wm3$w)
wm3dt = data.table(cat_wm3)
wm3dt[ , swrevf:=NA_real_]
wm3dt[ intersect(order(year), which(!is.na(w))), swrevf := cumprod(w), by=cowcode]
cat_wm3 <- as.data.frame(wm3dt)
summary(cat_wm3$swrevf)
sd(cat_wm3$swrevf)

############################################################
#### Outcome MSM with LHRPs (t+2), Weighting Model 3#######
############## Treatment=CmAT Review #####################
############################################################
msm3 <- lm(farissf2 ~ CATreview + CATnumprevrevs+ 
             CATyrssincereview + CATrevprevyr+ CATlintt, weights = swrevf, 
           data = cat_wm3)

######### Appendix Table A2 ######### 
stargazer (msm1, msm3)

###################################################################
###################################################################
### Appendix Table A3 #############################################
## Estimated CmAT review treatment and treatment history effects ##
##########Weighting Model 1 #######################################
###################################################################
### Data loading and pre-processing
##############################################
## Get the results directory
## Results for LHRPS (t+2), Weighting Model 1
model1<-read.csv("CAT_LHRPS__model_1_alpha_0_iterations_20000_outcomelead_2.csv")


model1_point <- model1[1,]
model1_95ci_lower <- apply(model1, 2, function(x) quantile(x, .025))
model1_95ci_upper <- apply(model1, 2, function(x) quantile(x, .975))
model1_50ci_mean <- apply(model1, 2, function(x) quantile(x, .50))


out1low<-t(model1_95ci_lower)
out1upper<-t(model1_95ci_upper)
g1<-rbind(model1_50ci_mean,out1low,out1upper)
ci1<-t(g1)
stargazer(ci1)

###################################################################
### Appendix Table A4 #############################################
## Estimated CmAT review treatment and treatment history effects ##
##########Weighting Model 3 #######################################
###################################################################
### Data loading and pre-processing
##############################################
model2<-read.csv("CAT_LHRPS__model_2_alpha_0_iterations_20000_outcomelead_2.csv")


model2_point <- model2[1,]
model2_95ci_lower <- apply(model2, 2, function(x) quantile(x, .025))
model2_95ci_upper <- apply(model2, 2, function(x) quantile(x, .975))
model2_50ci_mean <- apply(model2, 2, function(x) quantile(x, .50))

out2low<-t(model2_95ci_lower)
out2upper<-t(model2_95ci_upper)
g2<-rbind(model2_50ci_mean,out2low,out2upper)
ci2<-t(g2)
stargazer(ci2)

### Weighting Model 2 ###
# Denominator model (YEAR + treatment(review) history variables + time-varying covariates)
denom.p <- predict(wm2, type = "response")
# Numerator model (YEAR + treatment(review) history variables)
numer.fit <-  glm(CATreview ~  CATrevprevyr + CATyrssincereview + CATnumprevrevs + 
                    nevdem + yr, family="binomial", data=cat_wm2)
numer.p <- predict(numer.fit, type = "response")
## Calculate weights
## Calculate weights
d.pr.tr <- ifelse(cat_wm2$CATreview == 1, denom.p, 1-denom.p)
n.pr.tr <- ifelse(cat_wm2$CATreview == 1, numer.p, 1-numer.p)
cat_wm2$w <- n.pr.tr/d.pr.tr
summary(cat_wm2$w)
sd(cat_wm2$w)
wm2dt = data.table(cat_wm2)
wm2dt[ , swrevf:=NA_real_]
wm2dt[ intersect(order(year), which(!is.na(w))), swrevf := cumprod(w), by=cowcode]
cat_wm2 <- as.data.frame(wm2dt)
summary(cat_wm2$swrevf)
sd(cat_wm2$swrevf)

#### Outcome MSM with CIRI (t+2), Weighting Model 2#######
cat_wm2$tortf2 <- factor(cat_wm2$tortf2)
msm2 <- polr(tortf2 ~ CATreview + CATnumprevrevs+ CATyrssincereview+ CATrevprevyr+ 
               CATlintt ,method='probit',  weights = swrevf, 
             data = cat_wm2)

### Weighting Model 4 ###
# Denominator model (YEAR + treatment(review) history variables + time-varying covariates)
denom.p <- predict(wm4, type = "response")
# Numerator model (YEAR + treatment(review) history variables)
numer.fit <-  glm(CATreview ~  CATrevprevyr + CATyrssincereview + CATnumprevrevs + 
                    nevdem + yr, family="binomial", data=cat_wm4)
numer.p <- predict(numer.fit, type = "response")
## Calculate weights
## Calculate weights
d.pr.tr <- ifelse(cat_wm4$CATreview == 1, denom.p, 1-denom.p)
n.pr.tr <- ifelse(cat_wm4$CATreview == 1, numer.p, 1-numer.p)
cat_wm4$w <- n.pr.tr/d.pr.tr
summary(cat_wm4$w)
sd(cat_wm4$w)
wm4dt = data.table(cat_wm4)
wm4dt[ , swrevf:=NA_real_]
wm4dt[ intersect(order(year), which(!is.na(w))), swrevf := cumprod(w), by=cowcode]
cat_wm4 <- as.data.frame(wm4dt)
summary(cat_wm4$swrevf)
sd(cat_wm4$swrevf)

#### Outcome MSM with CIRI (t+2), Weighting Model 4 #######
cat_wm4$tortf2 <- factor(cat_wm4$tortf2)
msm4 <- polr(tortf2 ~ CATreview + CATnumprevrevs+ CATyrssincereview+ CATrevprevyr+ 
                CATlintt ,method='probit',  weights = swrevf, 
             data = cat_wm4)

####### ####### #######
####### Table A5  #######
####### ####### #######
stargazer(msm2, msm4)

####### ####### ############## ####### ############## ####### #######
####### Figure A5  ############## ####### ############## ####### #######
####### Simulated predicted probabilities, Weighting Model 2  #######
####### ####### ####### ####### ####### ############## ####### #######
## Load Models - base results

model1<-read.csv("CAT_CIRI__model_1_alpha_0_iterations_20000_outcomelead_2.csv")
m1_coeff<-rbind(model1$CATreview,model1$CATrevprevyr,model1$CATyrssincereview,model1$CATnumprevrevs,model1$CATlintt)

####Code to obtain the mean coefficients and their 95% CI
m1_95ci_upper <- apply(model1, 2, function(x) quantile(x, .975))
m1_95ci_lower <- apply(model1, 2, function(x) quantile(x, .025))
m1_mean<- apply(model1, 2, function(x) quantile(x, .50))
out1mean<-(m1_mean)
out1low<-(m1_95ci_lower)
out1upper<-(m1_95ci_upper)
g1<-rbind(out1mean,out1low,out1upper)
ci1<-t(g1)

xc.nrrev<-cbind(0, 0, 2.5, c(0:5), 13.5)
head(xc.nrrev)

ystar.nrrev<-(xc.nrrev)%*%(m1_coeff)
ystar.nrrev<-t(ystar.nrrev)
dim(ystar.nrrev)
m1_zeta<-cbind(model1$X0.1, model1$X1.2)
m1_zeta<-as.data.frame(m1_zeta)
###
pr1<- pnorm(m1_zeta$V1,mean=ystar.nrrev,sd=1)
pr2<- pnorm(m1_zeta$V2,mean=ystar.nrrev,sd=1)-pnorm(m1_zeta$V1,mean=ystar.nrrev,sd=1)
pr3<-1-pnorm(m1_zeta$V2,mean=ystar.nrrev,sd=1)

pr1mean<-apply(pr1,2,median)
pr2mean<-apply(pr2,2,median)
pr3mean<-apply(pr3,2,median)

postscript("FigureA5.eps", width=6.5, height=4, horizontal = FALSE, 
           onefile = FALSE)
par(cex.axis=1.5, cex.lab=1.2, cex.main=1.7, cex.sub=1.7)
plot(c(0:5),pr1mean,ylim=c(0,1),xlim=c(0,5),lwd=2,type="l",lty=3,
     xlab="Number of Previous Reviews ", ylab="Probability of observing a given level of Torture")

lines(c(0:5),pr2mean,ylim=c(0,1),lwd=2,lty=2,type="l")
lines(c(0:5),pr3mean,ylim=c(0,1),lwd=2,lty=4,type="l")
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(3,0,0),lwd=2, bty="n",cex=.9)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,2,0),lwd=2, bty="n",cex=.9)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,0,4),lwd=2, bty="n",cex=.9)

p1.low<- apply(pr1,2,quantile,prob=.025)
p1.high<- apply(pr1,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p1.low,  
         y1=p1.high,lwd=1) 
p2.low<- apply(pr2,2,quantile,prob=.025)
p2.high<- apply(pr2,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p2.low,  
         y1=p2.high,lwd=1) 
p3.low<- apply(pr3,2,quantile,prob=.025)
p3.high<- apply(pr3,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p3.low,  
         y1=p3.high,lwd=1)
dev.off()



####### ####### ############## ####### ############## ####### #######
####### Figure A6  ############## ####### ############## ####### #######
####### Simulated predicted probabilities, Weighting Model 4  #######
####### ####### ####### ####### ####### ############## ####### #######

model2<-read.csv("CAT_CIRI__model_2_alpha_0_iterations_20000_outcomelead_2.csv")
m2_coeff<-rbind(model2$CATreview,model2$CATrevprevyr,model2$CATyrssincereview,model2$CATnumprevrevs,model2$CATlintt)
m2_95ci_upper <- apply(model2, 2, function(x) quantile(x, .975))
m2_95ci_lower <- apply(model2, 2, function(x) quantile(x, .025))
m2_mean<- apply(model2, 2, function(x) quantile(x, .50))

out2mean<-(m2_mean)
out2low<-(m2_95ci_lower)
out2upper<-(m2_95ci_upper)
g2<-rbind(out2mean,out2low,out2upper)
ci2<-t(g2)


xc.nrrev<-cbind(0, 0, 2.63, c(0:5), 13.2)
head(xc.nrrev)

ystar.nrrev<-(xc.nrrev)%*%(m2_coeff)
ystar.nrrev<-t(ystar.nrrev)
dim(ystar.nrrev)
m2_zeta<-cbind(model2$X0.1, model2$X1.2)
m2_zeta<-as.data.frame(m2_zeta)
###
pr1<- pnorm(m2_zeta$V1,mean=ystar.nrrev,sd=1)
pr2<- pnorm(m2_zeta$V2,mean=ystar.nrrev,sd=1)-pnorm(m2_zeta$V1,mean=ystar.nrrev,sd=1)
pr3<-1-pnorm(m2_zeta$V2,mean=ystar.nrrev,sd=1)
pr1mean<-apply(pr1,2,median)
pr2mean<-apply(pr2,2,median)
pr3mean<-apply(pr3,2,median)



postscript("FigureA6.eps", width=6.5, height=4, horizontal = FALSE, 
           onefile = FALSE)
par(cex.axis=1.5, cex.lab=1.2, cex.main=1.7, cex.sub=1.7)
plot(c(0:5),pr1mean,ylim=c(0,1),xlim=c(0,5),lwd=2,type="l",lty=3,
     xlab="Number of Previous Reviews ", ylab="Probability of observing a given level of Torture")

lines(c(0:5),pr2mean,ylim=c(0,1),lwd=2,lty=2,type="l")
lines(c(0:5),pr3mean,ylim=c(0,1),lwd=2,lty=4,type="l")
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(3,0,0),lwd=2, bty="n",cex=.9)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,2,0),lwd=2, bty="n",cex=.9)
legend(1,1, c("No torture","Some torture","Widespread torture"),
       lty=c(0,0,4),lwd=2, bty="n",cex=.9)

p1.low<- apply(pr1,2,quantile,prob=.025)
p1.high<- apply(pr1,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p1.low,  
         y1=p1.high,lwd=1) 
p2.low<- apply(pr2,2,quantile,prob=.025)
p2.high<- apply(pr2,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p2.low,  
         y1=p2.high,lwd=1) 
p3.low<- apply(pr3,2,quantile,prob=.025)
p3.high<- apply(pr3,2,quantile,prob=.975)
segments(x0=seq(0,5,1), 
         y0=p3.low,  
         y1=p3.high,lwd=1)
dev.off()


####### ####### ############## ####### #######  #######
####### Table A6 - Standard Regression Models  #######
####### ####### ############## ####### #######  #######
rm(list=ls())
cat<- read.csv("CATReporting_Final.csv", header=TRUE)
cat <- subset(cat, cat$CATob==1)
cat$CATregiondensity <- 100*cat$CATregiondensity
cat$HRTpercent <- 100*cat$hrtprop
cat$yr <-as.factor(cat$year)
cat$lngdp <- log(cat$gdppc)
cat$lnpop <- log(cat$poptot)
catm1 <- subset(cat, !is.na(cat$gdppc) & !is.na(cat$poptot) & !is.na(cat$farissmean) & !is.na(cat$farissf2) &  
                !is.na(cat$polity) & !is.na(cat$CATreview)  & !is.na(cat$nhriest) 
              & !is.na(cat$CATArt22)  & !is.na(cat$CATregiondensity) & !is.na(cat$nevdem) &
                !is.na(cat$regional) & !is.na(cat$transdem) & 
                !is.na(cat$transition2) &  !is.na(cat$hrtprop))

catm2 <- subset(cat, !is.na(cat$gdppc) & !is.na(cat$poptot) & !is.na(cat$tort) & !is.na(cat$tortf2) &  
                  !is.na(cat$polity) & !is.na(cat$CATreview)  & !is.na(cat$nhriest) 
                & !is.na(cat$CATArt22)  & !is.na(cat$CATregiondensity) & !is.na(cat$nevdem) &
                  !is.na(cat$regional) & !is.na(cat$transdem) & 
                  !is.na(cat$transition2) &  !is.na(cat$hrtprop))
catm2$tortf2 <- as.factor(catm2$tortf2)

catm3 <- subset(cat, !is.na(cat$gdppc) & !is.na(cat$poptot) & !is.na(cat$farissmean) & !is.na(cat$farissf2) &  
                  !is.na(cat$polity) & !is.na(cat$CATreview)  & !is.na(cat$nhriest) 
                & !is.na(cat$CATArt22)  & !is.na(cat$CATregiondensity) & !is.na(cat$nevdem) &
                  !is.na(cat$regional) & !is.na(cat$transdem) & 
                  !is.na(cat$transition2) &  !is.na(cat$hrtprop) &  !is.na(cat$transparencyindex))

catm4 <- subset(cat, !is.na(cat$gdppc) & !is.na(cat$poptot) & !is.na(cat$tort) & !is.na(cat$tortf2) &  
                  !is.na(cat$polity) & !is.na(cat$CATreview)  & !is.na(cat$nhriest) 
                & !is.na(cat$CATArt22)  & !is.na(cat$CATregiondensity) & !is.na(cat$nevdem) &
                  !is.na(cat$regional) & !is.na(cat$transdem) & 
                  !is.na(cat$transition2) &  !is.na(cat$hrtprop) &  !is.na(cat$transparencyindex))
catm4$tortf2 <- as.factor(catm4$tortf2)


m1 <-zelig(farissf2 ~  lngdp + lnpop +nhriest +polity + nevdem +transdem+transition2 +  
             farissmean  + HRTpercent + CATArt22 + CATregiondensity+regional + 
             CATreview + CATnumprevrevs + CATrevprevyr+ CATyrssincereview +yr , model="normal", data=catm1)

m2 <- zelig(tortf2 ~  lngdp + lnpop +nhriest +polity + nevdem +transdem+transition2 +  
              tort  + HRTpercent + CATArt22 + CATregiondensity+regional + 
              CATreview + CATnumprevrevs + CATrevprevyr+ CATyrssincereview +yr , model="oprobit", data=catm2)

m3 <-zelig(farissf2 ~  lngdp + lnpop +nhriest + transparencyindex + polity + nevdem +transdem+transition2 +  
             farissmean  + HRTpercent + CATArt22 + CATregiondensity+regional + 
             CATreview + CATnumprevrevs + CATrevprevyr+ CATyrssincereview +yr , model="normal", data=catm3)

m4 <- zelig(tortf2 ~  lngdp + lnpop +nhriest + transparencyindex +polity + nevdem +transdem+transition2 +  
              tort  + HRTpercent + CATArt22 + CATregiondensity+regional + 
              CATreview + CATnumprevrevs + CATrevprevyr+ CATyrssincereview +yr , model="oprobit", data=catm4)

stargazer(m1, m2, m3, m4)

sink(NULL)

